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Summary 

The  Kim-Arakawa  orographic  gravity-wave  drag  parametrization  scheme,  which  is  a  component  of  the  US 
Navy’s  NOGAPS  ALPHA  (Navy  Operational  Global  Atmospheric  Prediction  System,  Advanced-Level  Physics 
and  High  Altitude),  is  extended  to  include  the  effects  of  orographic  anisotropy  and  low-level  flow  blocking.  The 
algorithms  to  calculate  the  orographic  statistics  needed  for  the  parametrization  are  also  revised.  The  extended 
scheme  is  evaluated  against  mountain  waves  explicitly  simulated  with  COAMPS®t  (Coupled  Ocean/Atmosphere 
Mesoscale  Prediction  System)  of  NRL  (Naval  Research  Laboratory). 

Mountain-wave  simulations  over  Boulder,  Colorado,  USA,  are  used  for  representing  realistic  situations 
of  different  wave  activity  including  severe  downslope  windstorms.  The  simulations  are  area-averaged  and 
interpolated  to  the  vertical  grid  of  NOGAPS,  and  are  used  as  the  input  to  the  extended  Kim-Arakawa  scheme. 
The  scheme  is  calibrated  by  comparing  the  parametrized  vertical  distribution  of  the  momentum  fluxes  with 
the  counterpart  obtained  from  the  explicit  mesoscale  simulations.  Overall,  the  calibrated  scheme  successfully 
represents  the  simulated  magnitudes  and  vertical  divergences  of  the  momentum  fluxes.  A  flow  regime  diagram  is 
constructed  utilizing  a  time  series  of  the  simulations  to  further  evaluate  the  parametrization.  The  robustness  of  the 
orographic  statistics,  together  with  an  approximate  method  to  improve  it,  are  also  addressed. 

KEYWORDS:  Blocked-layer  drag  Explicit  gravity-wave  simulation  Form  drag  Gravity-wave  drag 
Off-line  evaluation  Regime  diagram  Wave  breaking 


1.  INTRODUCTION 

The  effects  of  orography  can  be  represented  by  various  means  in  large-scale 
models  of  the  atmosphere,  and  their  inclusion  is  crucial  for  successful  simulation  and 
forecast  of  weather  and  climate.  Ever  since  the  first  generation  gravity-wave  drag 
(GWD)  parametrization  schemes  based  on  the  theories  of  two-dimensional  (2-D), 
linear,  hydrostatic,  stationary  mountain  waves  over  an  idealized  isolated  mountain  were 
introduced,  considerable  progress  has  been  made  in  improving  the  parametrizations. 
The  most  recent  examples  of  such  advancement  are  the  inclusion  of  the  effects  of 
orographic  anisotropy  and  low-level  flow  blocking. 

Until  recently,  the  effects  of  orography  for  the  stably  stratified  atmosphere  have 
been  Heated  by  separate  parametrizations:  the  ‘breaking’  of  unresolved  gravity  waves 
launched  by  subgrid-scale  orography  as  presented  by,  e.g.  Boer  et  al.  (1984),  Palmer 
et  al.  (1986)  and  McFarlane  (1987),  and  the  ‘blocking’  or  ‘stagnation’  of  low-level  flow 
through  enhanced  resolved  orography,  which  increases  planetary  wave  forcing  on  the 
large-scale  flow  (e.g.  Wallace  et  al.  1983;  Iwasaki  and  Sumi  1986;  Palmer  and  Mansfield 
1986;  Tibaldi  1986).  The  concept  of  ‘flow  blocking’  was  originally  introduced  into 
GWD  parametrization  in  the  mid  1980s.  For  a  large  (inverse)  Froude  number,  Frol, 
associated  with  a  relatively  high  mountain,  the  low-level  flow  is  blocked  by  the  mountain 
upstream,  and  the  effective  mountain  height  becomes  lower  than  its  actual  value  under 

*  Corresponding  author:  Naval  Research  Laboratory.  Marine  Meteorology  Division,  Stop  2,  7  Grace  Hopper 
Avenue,  Monterey,  CA  93943,  USA.  e-mail:  kimyj@nrlmry.navy.mil 
t  COAMPS  is  a  registered  trademark  of  the  Naval  Research  Laboratory. 

if  The  inverse  Froude  number  is  defined  as  Ftq  =  IiNq/Uo,  where  h  is  the  mountain  height,  and  No  and  Uq  are  the 
Brunt-Vaisala  frequency  and  the  horizontal  wind  magnitude,  respectively,  averaged  over  a  depth  which  defines 
low  levels.  For  GWD  parametrization,  h  is  usually  defined  as  a  multiple  of  the  standard  deviation  of  orography, 
cr/j ,  for  a  grid  box. 

©  Royal  Meteorological  Society,  2005. 
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the  2-D  framework  (see  Fig.  1,  reproduced  from  Kim  and  Arakawa  1994).  To  account 
for  this  effect,  Palmer  et  al.  (1986)  posed  an  upper  limit  on  the  height  (400  m),  while 
Pierrehumbert  (1986)  used  an  asymptotic  flux  function  that  gives  a  smooth  transition 
between  the  blocking  and  non-blocking  states;  see  Table  2  in  Kim  and  Arakawa  (1995, 
hereafter  KA95)  for  a  list  of  expressions  to  treat  2-D  flow  blocking.  This  concept  of 
limiting  the  height  of  orography  upstream  was  common  among  the  first-generation 
parametrization  schemes  before  the  3-D  nature  of  orography  was  fully  considered. 
In  other  words,  the  2-D  conceptualization  of  flow  blocking  considered  only  the  flow 
‘over’  the  mountain,  or  the  flow  being  blocked  completely  below  the  flow  separation  line 
that  occurs  on  the  mountain’s  flanks,  whereas  the  more  complete  3-D  conceptualization 
also  considers  also  the  flow  'around'  the  mountain  below  the  separation  line  (see  Fig.  1 
of  Lott  and  Miller  (1997),  hereafter  LM97). 

Surface  friction  is  known  in  general  to  reduce  mountain-wave  amplitudes  and  wave 
breaking  aloft  (e.g.  Richard  et  al.  1989;  Olafsson  and  Bougeault  1997,  hereafter  OB97; 
Doyle  et  al.  2000;  Leutbecher  and  Volkert  2000;  Doyle  and  Durran  2002;  Peng  and 
Thompson  2003).  Surface  friction  can  be  considered  more  important  than  GWD  with 
regard  to  momentum  exchange  with  the  surface  (e.g.  Shutts  and  Broad  1993)  and  is  often 
enhanced  to  represent  the  ‘form  drag’  through  increased  effective  surface  roughness 
due  to  turbulence  generated  by  subgrid-scale  orography  and  vegetation  under  neutral 
conditions  (e.g.  Wood  and  Mason  1993;  Milton  and  Wilson  1996,  hereafter  MW96; 
Belcher  and  Hunt  1998;  Gregory  et  al.  1998,  hereafter  GSM98).  This  form  drag  concept 
has  been  extended  to  a  scale  larger  than  turbulence,  where  it  was  originally  developed, 
to  represent  the  drag  generated  by  a  layer  of  blocked  flow  under  stable  conditions 
due  to  flow  past  subgrid-scale  orography  (e.g.  LM97;  Scinocca  and  McFarlane  2000, 
hereafter  SM00;  Webster  et  al.  2003;  Zadra  et  al.  2003;  Alpert  2004;  Brown  and  Webster 
2004).  As  also  pointed  out  by  GSM98,  the  traditional  definition  of  form  drag  refers  to 
turbulence-scale  orographic  roughness  under  neutral  or  unstable  conditions,  whereas 
drag  due  to  flow  blocking  is  to  account  for  subgrid-scale  (larger  than  turbulence-scale) 
blocking  or  splitting  of  flow  under  stable  conditions.  In  this  study,  therefore,  we  refer  to 
the  latter  process  as  ‘blocked-layer  drag’  (BLD).  As  a  more  physical  alternative  to  the 
resolved-scale  orographic  height  enhancement  discussed  above,  a  BLD  parametrization 
is  now  implemented  in  a  subgrid-scale  sense,  providing  an  additional  source  of  low-level 
drag  under  the  3-D  framework,  and  is  often  integrated  as  a  part  of  GWD  parametrization 
such  as  by  LM97  and  SM00. 

The  scheme  proposed  by  KA95  parametrizes  the  drag  due  to  breaking  and  trapping 
of  both  hydrostatic  and  non-hydrostatic  gravity  waves,  and  distinguishes  between  the 
2-D  flow-blocking  ‘upstream’  state  and  wave  breaking  ‘downstream’  state  based  on 
Pierrehumbert’s  (1986)  formulations,  without  incurring  drag  enhancement  due  to  flow 
blocking  in  the  3-D  sense.  That  is,  the  blocking  due  to  stagnant  flow  formed  upstream 
of  orography  reduces  the  effective  height  of  orography,  resulting  in  the  decrease  of 
drag,  whereas  low-level  wave  breaking  (LLWB)  or  lee-wave  trapping  in  the  downstream 
region  involves  a  significant  increase  in  drag  when  certain  flow  conditions  are  met 
(Fig.  1).  This  scheme  collectively  treats  any  enhancement  of  GWD  at  low  levels  due 
to  LLWB  or  lee-wave  trapping  in  terms  of  the  2-D  resonant  amplification  of  GWD 
(Peltier  and  Clark  1979);  it  is  physically  similar  to  that  of  GSM98,  but  considers 
only  partial  effects  of  orographic  anisotropy  through  additional  orographic  statistics. 
Its  implementations  into  large-scale  models  are  reported  in  Kim  (1996;  hereafter  K96), 
Alpert  et  al.  (1996)  and  also  Kim  and  Hogan  (2004). 

The  primary  goal  of  the  present  study  is  to  extend  and  evaluate  the  KA95  scheme 
by  including  3-D  effects  of  orography  and  also  by  including  a  BLD  formulation. 
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Figure  1.  The  key  processes  that  the  orographic  gravity- wave  drag  (GWD)  scheme  of  Kim  and  Arakawa  (1995) 
attempts  to  parametrize.  In  the  downstream  region,  low-level  wave  breaking  and/or  wave  trapping  of  lee  waves  can 
contribute  to  GWD.  The  vertical  dotted  line  is  described  in  appendix  A.  (Taken  from  Kim  and  Arakawa  (1994).) 


We  introduce  a  new  orographic  statistical  parameter  to  consider  the  effect  of  oro¬ 
graphic  anisotropy.  We  derive  the  BLD  formulation  basically  following  earlier  studies 
by  LM97  and  SMOO,  which  are  based  on  numerical  studies  of  3-D  flow  past  isolated 
mountains  (e.g.  Hunt  and  Snyder  1980;  Smolarkiewicz  and  Rotunno  1989;  Miranda  and 
James  1992,  hereafter  MJ92;  Schar  and  Smith  1993;  Olafsson  and  Bougeault  1996, 
hereafter  OB96;  OB97;  and  Schar  and  Durian  1997).  A  rigorous  evaluation  of  the 
parametrization  with  direct  measurements  is  usually  difficult  due  to  limited  availability 
and  quality  of  high-resolution  observation  data.  An  alternative  approach  to  evaluating 
the  parametrization  with  observations  is  the  use  of  a  high-resolution  mesoscale  model, 
which  has  been  evaluated  by  various  means,  e.g.  against  observational  data  and/or  other 
verified  model  simulations  as  by,  e.g.  KA95,  Broad  (1996),  LM97  and  GSM98.  In  this 
study  we  evaluate  the  extended  KA95  parametrization,  hereafter  KA95+,  using  explicit 
simulations  of  dry  mountain  waves  obtained  from  a  3-D  mesoscale  model. 

In  section  2,  we  present  the  reformulated  versions  of  orographic  statistics  origi¬ 
nally  developed  for  the  KA95  scheme  and  also  a  new  parameter  for  taking  into  account 
orographic  anisotropy.  In  section  3,  we  summarize  the  KA95+  scheme;  section  4 
presents  an  evaluation  of  the  scheme.  We  first  describe  our  mesoscale  model  that  explic¬ 
itly  simulates  mountain  waves,  and  present  case-studies  of  the  simulations  representing 
typical  mountain  waves.  We  then  evaluate  the  scheme  using  the  cases  by  comparing 
simulated  and  parametrized  results.  We  construct  a  regime  diagram  using  all  times  of 
the  simulations  in  an  effort  to  further  evaluate  the  scheme.  We  also  discuss  important 
issues  regarding  the  robustness  of  the  orographic  statistics  and  the  effects  of  moisture. 
Concluding  remarks  on  remaining  issues  of  the  KA95+  scheme  in  particular  and  those 
of  the  GWD  parametrization  in  general  are  given  in  section  5.  Appendix  A  describes  a 
non-local  version  of  the  parametrization  scheme  and  appendix  B  introduces  a  method  to 
obtain  weighted  averages  of  the  orographic  statistics. 


2.  Statistics  of  subgrid-scale  orography  for  the  KA95+  scheme 

The  KA95  scheme  requires  higher  order  statistics  of  orography  as  well  as  the  stan¬ 
dard  deviation  (07,)  of  subgrid-scale  orographic  heights.  It  utilizes  the  semi-empirical  but 
geophysical  relationships  between  the  configuration  of  subgrid-scale  orography  and  the 
physical  characteristics  of  the  corresponding  flow,  which  were  obtained  from  extensive 
2-D  mesoscale  mountain-wave  simulations  involving  over  100  cases  with  various  shapes 
and  sizes  of  mountains.  Here,  we  discuss  these  orographic  statistics. 
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(a)  Orographic  asymmetry 

Lee  waves  that  are  trapped  near  the  surface  may  propagate  laterally  to  neighbouring 
grid  cells,  and  similarly  vertically  propagating  waves  may  traverse  out  of  the  origi¬ 
nal  model  column  to  neighbouring  columns.  Thus,  GWD  originating  from  a  grid  cell 
with  orography  may  even  influence  a  neighbouring  cell  with  flat  topography.  As  the 
resolution  of  models  increases,  this  non-local  nature  of  drag  becomes  more  of  a  criti¬ 
cal  issue;  however,  current  GWD  parametrizations  that  assume  a  local  column  physics 
framework  ignore  this  effect.  This  issue  was  addressed  in  Kim  (1992)  by  introducing 
a  new  parameter  which  considers  the  effect  of  neighbouring  grid  domains  (appendix  A 
describes  this  issue  further).  This  approach,  however,  introduces  significant  computa¬ 
tional  inefficiency  under  current  parallel  computing  architecture  due  to  the  need  for 
communication  with  neighbouring  grid  cells.  Therefore,  as  a  first  step  toward  addressing 
the  larger  issue  of  non-local  drag  variance,  KA95  incorporated  this  effect  within  a  grid 
box,  which  is  revisited  here  (for  more  discussion  see  appendix  A). 

The  ‘orographic  asymmetry’  (CM)  measures  the  asymmetry  of  subgrid-scale  orog¬ 
raphy  and  its  relative  location  to  the  model  grid  box  (see  appendix  B  of  KA95  for  the 
original  definition).  We  generalize  the  algorithm  to  calculate  OA  for  any  orographic  data 
of  any  resolution,  with  either  fixed  intervals  for  grid  point  models  or  variable  intervals 
for  the  Gaussian  grid  of  spectral  models: 

CM  =  1  —  (Nd/Nu),  (1)* 

where  Nu  and  ND  denote,  respectively,  the  numbers  of  subgrid-scale  high-resolution 
orographic  elevations  in  the  upstream  and  downstream  regions  (relative  to  the  wind 
direction)  higher  than  the  grid-box  average  of  the  orographic  heights  (which  is  a 
good  approximation  for  high  resolutions  to  the  ‘mode’  used  in  K96,  Eq.  (B.l)).  The 
calculations  of  Au  and  N D  for  the  four  representative  wind  directions  (west,  south, 
south-west  and  north-west)  arc  illustrated  in  Fig.  2.  Depending  on  the  pre-determined 
directions  of  the  low-level  wind  (zonal,  meridional,  or  diagonal  in  either  direction),  the 
regions  of  the  upstream  and  downstream  areas  arc  defined  as  shown  in  the  figure.  OA  is 
obtained  using  Eq.  (1)  for  the  four  representative  directions,  and  OA  (east,  north,  north¬ 
east  and  south-east)  =  —CM  (west,  south,  south-west  and  north-west),  respectively  (e.g. 
CM  (east)  is  identical  to  —CM  (west)).  In  a  large-scale  model,  CM  is  selected  at  each  time 
step  depending  on  the  wind  direction. 

If  the  mountain  is  symmetric  and  located  in  the  centre  of  the  grid  domain,  CM  is 
zero.  Although  the  mountain  is  located  in  the  centre  of  domain,  however,  CM  becomes 
positive  (negative)  if  the  mountain  is  skewed  to  downstream  (upstream)  as  illustrated  in 
Figs.  3(a)  and  (b).  This  design  is  based  on  the  2-D  mountain-wave  simulations  by  KA95 
and  others  showing  that  steeper  slopes  on  the  leeward  sides  involve  stronger  nonlinear 
waves  and  more  likelihood  for  LLWB  and/or  non-hydrostatic  wave  trapping.  CM  can 
also  distinguish  between  the  two  configurations  of  orography  which  involve  the  same 
mountain  in  the  grid,  but  different  locations  relative  to  the  grid  (see  Figs.  3(c)  and  (d)) — 
one  case  including  mostly  downstream  of  the  mountain  (CM  >  0)  and  the  other  upstream 


*  We  limit  the  range  of  OA  as  in  KA95  to  prevent  peculiarly  large  values  due  to  a  small  number  of  orographic  data 
points  in  a  grid  cell,  so  that  —  1  OA  ^  1,  which  is  typically  found  with  real  orographic  data  of  any  resolution. 
An  alternative  formula  has  been  suggested  (Dr  N.  Wood,  personal  communication),  which  does  not  require  the 
arbitrary  limit  and  ensures  the  anti-symmetry  of  OA  for  opposite  wind  directions  even  with  a  small  number  of 
orographic  data  points:  OA  =  c(Nv  —  ND)/(NV  +  ND),  where  the  case  with  c  =  3  is  then  consistent  with  the 
original  definition  of  Eq.  (1). 
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Figure  2.  Method  of  dividing  a  grid  box  for  calculating  the  numbers  of  the  subgrid-scale  orographic  height 
data  points  in  the  upstream  (IV  ,  shaded)  and  downstream  (ND)  regions  to  be  used  for  obtaining  the  orographic 
asymmetry  ( OA ),  depending  on  the  low-level  wind  direction,  i.e.  (a)  westerly,  (b)  southerly,  (c)  south-westerly, 
and  (d)  north-westerly.  For  the  diagonal  wind  directions  ((c)  and  (d)),  two  side  subgrid  boxes  are  summed  and 
halved  to  be  used  both  for  the  upstream  and  downstream  regions.  (Compare  with  Fig.  6.) 


(a)  OA  >  0 


0  Ax/2  Ax  0  Ax/2  Ax 

(c)  OA  >  0 


0  Ax/2  Ax  0  Ax/2  Ax 

Figure  3.  Illustration  of  orographic  asymmetry  (OA)  for  an  isolated  mountain  sharing  the  same  value  of  cr/, 
(standard  deviation  of  topography)  in  the  grid  domain  of  side-length  Ax.  For  a  mountain  that  is  symmetric  and  in 
the  centre  (not  shown),  OA  is  defined  as  zero.  Although  the  mountain  is  centred,  if  it  is  skewed  to  (a)  downstream 
((b)  upstream)  OA  is  (a)  positive  ((b)  negative).  If  the  mountain  is  symmetric,  but  located  toward  the  (c)  upwind 
((d)  downwind)  direction  OA  is  (c)  positive  ((d)  negative).  Note  that  2-D  cases  are  shown  here  for  simplified 
visualization,  but  the  actual  statistics  are  calculated  for  the  grid-box  area  rather  than  for  grid  length  Ax. 


(d)  OA  <  0 


(b)  OA  <  0 


(OA  <  0).  This  is  also  based  on  results  of  the  simulations  by  KA95  and  others  show¬ 
ing  that  downstream  regions  are  more  likely  to  contain  stronger  LLWB  and/or  non¬ 
hydrostatic  wave  trapping.  The  two  idealized  configurations  shown  in  Figs.  3(c)  and  (d) 
share  the  same  mountain  in  the  same  domain  and  thus  are  distinguishable  not  by  the 
value  of  <7h  in  the  grid  but  by  the  values  of  OA  with  opposite  sign.  In  the  parametrization, 
a  larger  OA  generally  corresponds  to  a  stronger  GWD. 
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(a)  OC  -  Larger  (b)  OC  -  Smaller 


0  Ax/2  Ax  0  Ax/2  Ax 


Figure  4.  Illustration  of  orographic  convexity  (OC)  for  an  idealized,  symmetric  and  isolated  mountain  case  in 
the  grid  domain  of  side-length  Ax.  A  (a)  sharper  ((b)  duller)  mountain  corresponds  to  (a)  larger  ((b)  smaller) 
OC.  Note  that  2-D  cases  are  shown  here  for  simplified  visualization,  but  the  actual  statistics  are  calculated  for  the 
grid-box  area  rather  than  for  the  grid  length  (Ax). 


( b )  Orographic  convexity 

The  ‘orographic  convexity’  (OC)  represents  the  sharpness  (and  slope)  of  the  moun¬ 
tain^),  which  is  linked  to  the  characteristics  of  the  corresponding  mountain  wave.  The 
expression  for  OC  taken  from  KA95  is: 

Nx 

OC=(l/Nxa£)J2(hi-h)\  (2) 

i=l 

where  x  represents  the  horizontal  direction,  Nx  denotes  the  number  of  subgrid-scale 
orographic  height  points  within  the  large-scale  grid,  and  the  overbar  denotes  the  large- 
scale  grid  average.  The  sharpness  of  mountain  has  been  considered  earlier  in  the 
parametrization  by  Phillips  (1984).  In  general,  a  sharper  (duller)  mountain  is  associated 
with  larger-  (smaller-)  amplitude  waves  (Fig.  4).  As  shown  by  KA95,  OC  is  well 
correlated  with  the  mountain’s  ‘vertical  aspect  ratio’  (height  to  half-width;  see  Fig.  16 
of  KA95)  in  the  case  of  isolated,  symmetric  mountains,  which  characterizes  whether 
the  waves  are  launched  with  characteristics  of  hydrostatic  (low  vertical  aspect  ratio)  or 
non-hydrostatic  (high  vertical  aspect  ratio)  waves. 

(c)  Non-dimensional  effective  orographic  length 

The  ‘effective  orographic  length’  ( Lx )  is  the  subgrid-scale  mountain  width  mea¬ 
sured  at  the  critical  orographic  height  (hc),  which  is  integrated  over  the  grid  and  nor¬ 
malized  by  the  grid  size  (see  Fig.  5).  Lx  was  designed  to  complement  Fro,  which  alone 
cannot  accurately  measure  the  nonlinearity  of  the  flow  over  complex  mountains  (KA95), 
by  estimating  the  bulk  volume  of  subgrid-scale  orography  that  is  associated  with  non¬ 
linearity  of  the  flow.  Strongly  nonlinear  flow  with  larger  Fro  involves  smaller  hc,  and 
thus  larger  Lx.  K96  used  the  expression,  hc  =  I  1 16.2  —  O.878o>,,  which  was  empirically 
derived  from  the  2-D  mountain-wave  simulations  of  KA95,  to  avoid  calculating  hc  in 
every  time  step  using  the  original  definition,  rr/,Frc/Fro,  where  Frc(~ 0.8)  is  a  prescribed 
critical  Froude  number.  In  this  study  we  further  assume,  based  on  some  earlier  tests,  that 
the  grid-box  average  of  the  subgrid-scale  orographic  heights  is  a  crude  approximation 
to  the  empirical  value  of  hc  and  now  use  the  average  in  calculating  Lx. 

The  expression  for  Lx  described  by  KA95  is  generalized  here  as: 

Lx  =  J2  L*</Ax  =  Nw/Nt.  (3) 

i 

Here,  Av  is  the  horizontal  grid  size,  Nw  denotes  the  number  of  subgrid-scale  orographic 
heights  (high-resolution  elevation  data  points)  along  the  centre  area  with  respect  to 
the  wind  crossing  over  four  one-eighth  grid  boxes  (for  the  directions  west  and  south; 
Figs.  6(a)  and  (b))  or  two  one-quarter  grid  boxes  (for  south-west  and  north-west; 
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Figure  5.  Illustration  of  the  non-dimensional  effective  orographic  length  (Lx  =  JA  LXJ Ax),  which  is  the  sum 
of  the  non-dimensionalized  horizontal  lengths  intersecting  the  subgrid-scale  mountain  at  the  critical  height  (hc) 
in  the  grid  of  length  Ax.  hc  is  approximated  in  this  study  by  the  grid-box  (area)  averaged  value  of  subgrid-scale 
orographic  heights.  The  value  of  Lx  is  calculated  for  the  four  representative  low-level  wind  directions  (described 
in  Fig.  6).  Note  that  here  2-D  cases  are  shown  for  simplified  visualization,  but  the  actual  statistics  are  calculated 

for  the  grid-box  area  rather  than  for  (Ax). 


Figure  6.  Methods  of  dividing  a  grid  box  for  calculating  the  number  of  subgrid-scale  orographic  height  data 
points  in  the  grid  box  higher  than  the  grid  average  ((Vw)  to  be  used  for  obtaining  the  effective  orographic  length, 
Lx.  The  total  number  of  orographic  height  data  points  ((VT),  which  is  also  needed  to  obtain  Lx  (see  Eq.  (3)),  is 
counted  for  the  same  subgrid-scale  boxes  as  for  (Vw.  Depending  on  the  low-level  wind  direction,  four  one-eighth 
boxes  or  two  one-quarter  boxes  (shaded)  are  regarded  as  the  centre  section  of  the  box  in  the  direction  of  the  flow. 


Figs.  6(c)  and  (d)),  which  arc  higher  than  the  grid  average  (both  considering  effectively 
the  half  of  the  grid  box  area  as  shaded  in  the  figure).  NT  denotes  the  total  number 
of  orographic  data  points  in  the  grid  box.  The  calculation  is  performed  for  the  four 
representative  wind  directions  similarly  to  that  of  OA,  and  Lx  (east,  north,  north-east 
and  south-east)  =  Lx  (west,  south,  south-west  and  north-west),  respectively.  In  a  large- 
scale  model,  Lx  is  selected  at  each  time  step  depending  on  the  wind  direction. 
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Figure  7.  Illustration  of  the  length  of  the  large-scale  grid  (A*)  and  the  subgrid-scale  effective  orographic  length 
(Lx)  in  the  direction  of  the  low-level  wind  (t/o)  for  an  isolated  mountain;  also  the  lengths  in  the  direction 
perpendicular  to  the  low-level  wind  (A^  and  L)r,  respectively).  The  critical  orographic  height,  hc,  which  is 
approximated  by  the  grid  mean  orographic  height  in  this  study,  is  shown  by  the  dashed  line.  The  orographic 
direction  (OD),  defined  by  Eq.  (4),  is  calculated  using  Lx  and  Lx h 


( d )  Orographic  direction 

It  is  known  that  3-D  spreading  of  mountain  waves  tends  to  reduce  wave  momentum 
fluxes  and  GWD  above  (e.g.  Shutts  1998).  It  was  reported  that  idealized  3-D  orography 
generates  only  about  the  half  of  the  momentum  flux  in  comparison  with  correspond¬ 
ing  2-D  orography  (Nappo  and  Chimonas  1992).  Orographic  anisotropy  was  earlier 
included  in  the  GWD  parametrizations  by  Baines  and  Palmer  (1990)  and  Shutts  (1990). 
The  KA95  scheme  was  originally  calibrated  based  on  a  2-D  framework  using  simula¬ 
tions  with  a  2-D  mountain-wave  model.  The  3-D  effects  of  orography  were  considered 
only  partially  through  OA  and  Lx,  which  were  calculated  for  the  eight  representative 
wind  directions  as  described  above.  Thus  the  KA95  scheme  did  not  fully  consider  the 
reduction  of  wave  amplitude  due  to  orographic  anisotropy.  As  a  result,  the  low-level 
Froude  number  (Fro  =  hNo/Uo),  which  is  one  of  the  major  measures  of  nonlinearity  of 
the  flow  for  the  parametrization,  inherently  assumed  the  mountain  of  infinite  length  in 
the  cross-wind  direction.  To  rectify  this  deficiency,  we  introduce  ‘orographic  direction’ 
(OD)  representing  the  orographic  anisotropy: 


(4) 


where  superscript  ±  denotes  the  cross-wind  direction,  i.e.  Ljr  denotes  Lx  for  the  direc¬ 
tion  perpendicular  to  the  low-level  wind  (see  Fig.  7).  OD  is  equivalent  to  the  mountain’s 
‘horizontal  aspect  ratio’  (cross-width  to  along-width)  or  inverse  ‘eccentricity’  (SM00) 
for  a  single  symmetric  mountain,  but  is  defined  more  generally  here  as  dominant  bulk 
subgrid-scale  orography  summed  over  the  grid  box  with  respect  to  the  representative 
wind  directions.  The  Froude  number  is  accordingly  redefined  as: 


No 

Fro  =  h — -OD, 
U0 


where  the  orographic  height,  h,  is  defined  as  2 07,  in  this  study. 


(5) 


(e)  Orographic  statistics  database 

The  orographic  statistics  systematically  consider  the  details  of  subgrid-scale  oro¬ 
graphy  (e.g.  shape,  size,  number,  distribution,  direction,  etc).  These  statistics  are 
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designed  to  be  calculated  offline  for  computational  efficiency.  First,  OC  and  07,,  which 
do  not  depend  on  the  wind  direction,  are  calculated.  Then,  lookup  tables  of  CM,  OD 
and  Lx  with  respect  to  the  representative  four  low-level  wind  directions  (west,  south, 
south-west  and  north-west)  arc  constructed.  For  consistency,  these  orographic  statistics 
derived  from  the  revised  algorithms  have  been  compared  for  the  2-D  ridge  cases  with 
those  from  the  original  algorithms  of  KA95.  When  used  online  in  a  large-scale  model, 
the  parametrization  scheme  first  determines  the  direction  of  the  low-level  wind  at  each 
time  step,  and  then  simply  looks  up  the  pre-calculated  tables  to  find  the  corresponding 
orographic  statistics. 


3.  Extension  of  the  KA95  orographic  GWD  parametrization  scheme 


(a)  The  gravity-wave  drag  parametrization 
The  formulations  for  the  GWD  parametrization  are  formally  the  same  as  those 
of  KA95  and  K96  except  that  Ftp  is  now  multiplied  by  OD  to  include  the  effects  of 
orographic  anisotropy  as  shown  in  Eq.  (5).  The  GWD  (r)  at  the  reference  level  (/iref)  is: 


m  l£/()|3 
tgwd  =  Po  E  - — G— — , 

Aeff  Aq 


E  =  (OA  +  2) 


CEFr0/Frc 


m  =  (  1  +  Lx) 


04+1 


G  = 


Ft2 

r'o 


Ft2  +  CGOC~l  ’ 


(6) 

(7) 


where  p  is  the  density;  N  is  the  Brunt-Vaisala  frequency,  U  is  the  horizontal  wind  speed 
projected  to  the  direction  of  the  low-level  wind;  E  is  the  ‘enhancement  factor’  applied 
only  at  the  reference  level  to  represent  the  nonlinear  enhancement  of  drag  by  resonant 
amplification  (e.g.  Peltier  and  Clark  1979)  in  the  downstream  regions  due  to  LLWB 
and/or  wave  trapping,  and  is  controlled  by  the  shape  and  location  of  subgrid-scale 
orography  in  the  grid  (OA)  and  the  nonlinearity  of  the  grid-scale  flow  (Ftp)  normalized 
by  its  critical  value;  m  is  the  ‘number  of  mountains’,  which  estimates  the  bulk  volume  of 
subgrid-scale  orography  associated  with  the  nonlinearity  of  the  flow  (Lx),  and  depends 
also  on  the  shape  and  location  of  subgrid-scale  mountain(s)  relative  to  the  grid  (OA);  G 
is  an  asymptotic  function  that  provides  a  smooth  transition  between  2-D  non-blocking 
and  blocking  cases  as  used  by  Pierrehumbert  (1986)  and  includes  the  influence  of  the 
vertical  mountain  aspect  ratio  through  OC,  empirically  applying  the  original  idea  by 
Pierrehumbert  (1986;  Eq.  (3.8));  and  /.eff  is  the  effective  grid  length,  which  was  set  to 
the  length  of  the  grid  box  in  K96,  but  can  be  used  practically  as  a  tuning  coefficient. 
We  set  Ce  =  0.8  and  Cq  =  0.5  as  originally  calibrated  with  mesoscale  simulations 
(KA95).  The  subscript  o  denotes  a  low-level  average,  which  in  this  study  is  between 
the  surface  and  2 07,  (=  h)  differing  from  the  original  definition  in  K96  of  the  depth 
of  the  atmospheric  boundary  layer.  The  expressions  in  Eq.  (7)  were  derived  semi- 
empirically  based  on  physical  concepts  and  evaluated  against  extensive  2-D  mountain- 
wave  simulations  (KA95). 

A  nonlinear  extension  of  the  original  Pierrehumbert  formulation,  in  ■  E,  becomes 
larger  as  OA  increases,  in  order  to  be  consistent  with  the  simulations  of  KA95  that 
indicate  downstream  regions  arc  more  likely  to  include  stronger  wave  activity,  and 
thus  a  greater  chance  for  drag  enhancement.  KA95  demonstrated  that  m  ■  E  accurately 
represents  at  low  levels  (see  Figs.  18(c)  and  25  of  KA95)  the  vertical  gradient  of  the 
Scorer  parameter  (Scorer  1949),  which  can  be  approximated  as  t2~N2/U2  (more 
rigorous  definitions  are  given  by,  e.g.  KA95  footnote  7,  or  Nance  and  Durran  (1998) 
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Eq.  (4)).  Analytical  and  numerical  simulation  studies  of  mountain  waves  show  that  the 
Scorer  parameter  influences  the  linear  resonant  flapping  of  non-hydrostatic  waves  (lee 
waves)  at  low  levels  (e.g.  Bretherton  1969;  Shutts  and  Broad  1993;  Doyle  et  al.  2002). 
KA95  further  found  that  it  also  serves  as  a  good  measure  of  nonlinear  resonant  LLWB. 
Whether  of  linear  or  nonlinear  nature,  this  wave  activity  at  low  levels  contributes  to 
increased  vertical  divergence  of  GWD.  KA95  demonstrated  that  the  vertical  gradient 
of  the  Scorer  parameter  cannot  be  uniquely  parametrized  by  Fro — the  main  flow 
parameter  used  in  earlier  conventional  GWD  schemes — for  different  situations,  say, 
with  or  without  LLWB  (see  Lig.  13  of  KA95).  Peng  and  Thompson  (2003)  also 
questioned  whether  Fro  in  GWD  parametrizations  could  properly  represent  boundary- 
layer  processes. 

A  model  layer  is  considered  unstable  if  the  minimum  Richardson  number  defined 


by: 


m  -  Frd ) 

(1  +  yfRl  Fro)2 


(8) 


becomes  less  than  its  critical  number,  i.e.  if  Rim  <  Ric  (typically,  Ric  —  0.25),  where  Ri 
is  the  (mean)  Richardson  number,  (defined  by  Eq.  (2.11)  in  KA95)  and  Fro  =  hoN/U. 
Here,  the  vertical  displacement  height,  ho,  is  expressed  combining  the  expressions  by 
Palmer  et  al.  (1986)  and  Pierrehumbert  (1986)  as: 


(h2o)i  = 


Ax  r,_|_  i 
m  PiNjUj ’ 


(9) 


where  subscript  i  denotes  the  vertical-layer  index  decreasing  upward.  The  expression 
for  the  critical  value  of  ho  is  then  derived  by  substituting  Ric  =  0.25  into  Eq.  (8)  as: 


(h  d)c  = 


U 

N 


(10) 


The  profile  of  GWD  is  determined  according  to  the  following  steps  (see  KA95  for 
more  details): 


•  The  reference -level  drag,  tgwd.  is  calculated  from  Eq.  (6). 

•  The  vertical  displacement  height,  ho,  is  computed  by  Eq.  (9). 

•  The  layer  stability  is  checked  by  evaluating  Rim  using  the  computed  ho  in  Eq.  (8). 

•  If  Rim  >  Ric,  t  is  unchanged  for  the  next  model  layer,  whereas,  if  Rim  si  Ric, 
z  is  calculated  from  Eq.  (9)  (with  i  +  1  replaced  by  /)  using  (ho)c  obtained  from 
Eq.  (10). 

It  was  shown  by  KA95,  based  on  2-D  mountain-wave  simulations,  that  LLWB 
or  wave  trapping  is  not  properly  represented  by  the  saturation  hypothesis  of  Lindzen 
(1981),  but  is  better  represented  by  the  vertical  gradient  of  the  Scorer  parameter.  For 
this  reason,  KA95  determined  the  vertical  decrease  of  the  drag  at  low  levels  (below  the 
‘interface  height’,  hmt)  by  the  ratio  of  the  Scorer  parameter  instead  of  the  saturation 
hypothesis  when  the  Ri  criterion  is  met  in  a  model  layer  as: 

—5—  =  Min  (  Ci  ,  1 

where  Q  =  1,  and  hmt  is  now  defined  as  the  level  where  Ri  first  decreases  with  height 
above  the  reference  level.  Moreover,  this  ratio  is  now  applied  regardless  of  04,  unlike 
K96  who  used  the  ratio  only  for  OA  >  0,  in  order  to  ensure  a  smoother  transition 
between  the  upstream  and  downstream  configurations. 
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(b)  The  blocked-layer  drag  parametrization 
In  the  introduction,  we  discussed  the  recent  advances  in  the  GWD  parametrization 
community  that  includes  an  extra  source  of  drag  due  to  flow  blocking  at  low  levels  in 
place  of  enhanced  resolved  orography.  Along  these  lines,  we  extend  KA95  by  including 
BLD  in  the  drag  parametrization.  Our  formulation  basically  follows  previous  pioneering 
parametrization  studies  (e.g.  LM97,  SM00)  and  takes  the  following  bulk  aerodynamic 
drag  form,  which  is  based  in  paid  on  scale  analysis  and  is  particularly  similar  to  that  of 
SM00;  however,  some  key  parameters  are  calculated  using  distinctly  different  methods: 

^bld  =  -po-£2CdA^L^hB\Uo\2,  (12) 

where  is  the  grid  box  area,  Cd  is  a  bulk  drag  coefficient  of  order  unity  (we  set  Cd  —  1 
in  this  study),  is  the  length  of  large-scale  grid  in  the  cross-wind  direction,  and  L^r  is 
as  used  in  Eq.  (4)  the  width  of  dominant  subgrid-scale  orography  along  the  cross-wind 
direction,  which  is  approximated  as  the  width  of  orography  measured  at  the  critical 
orographic  height  (see  Fig.  7);  h b  is  the  height  of  blocked  layer  defined  as: 

h  B  =  (Fro  ~  Frc)  >0.  (13) 

N0 

The  orographic  anisotropy  is  considered  through  /? b  that  includes  Fro  (Eq-  (13)), 
which  itself  is  multiplied  by  OD  (Eq.  (5)).  The  BLD  is  calculated  using  Eq.  (12) 
when  Fro  >  Frc,  and  applied  to  the  lowest  model  level  above  the  surface  and  linearly 
decreased  in  height  as  in  other  studies.  (More  details  are  given  in  subsection  4(c).) 


4.  Evaluation  of  the  KA95+  scheme  with  explicit  mountain-wave 

SIMULATIONS 

(a)  Mesoscale  model  COAMPS® 

For  explicit  simulations  of  mountain  waves,  we  use  the  atmospheric  component  of 
COAMPS®  (Hodur  1997).  This  mesoscale  model  uses  finite-difference  approximations 
to  describe  the  fully  compressible,  non-hydrostatic  equations  that  govern  atmospheric 
motions,  and  a  terrain-following  vertical  coordinate.  The  prognostic  variables  of  the 
model  are  the  horizontal  velocity,  vertical  velocity,  perturbation  Exner  function,  po¬ 
tential  temperature,  turbulent  kinetic  energy,  and  the  concentrations  of  moisture,  cloud 
drops,  ice  crystals,  snow,  rain,  graupel  and  aerosols.  The  horizontal  and  vertical  advec- 
tion,  pressure  gradient,  and  divergence  are  represented  by  second  order  accurate  dif¬ 
ferencing  in  this  application.  Nonlinear  instability  is  suppressed  by  a  fourth-derivative 
hyper-diffusion.  Turbulent  processes  are  represented  by  a  level-2.5  closure  parametriza¬ 
tion,  similar  to  the  level- 1.5  turbulence  closure  scheme  used  in  KA95’s  dry  mountain- 
wave  model,  which  was  essential  for  proper  development  of  mountain  waves  with  large 
amplitudes.  Lateral  boundary  conditions  make  use  of  the  Navy  Operational  Global 
Atmospheric  Prediction  System  (NOGAPS,  Hogan  and  Rosmond  1991)  forecast  fields 
following  Davies  (1976).  The  nested-grid  boundary  conditions  are  formulated  using  a 
one-way-interaction  approach  with  a  horizontal  resolution  of  each  nested  mesh  that  is 
1  /3  of  the  parent  grid  mesh. 

Shutts  and  Broad  (1993)  and  Broad  (1996)  performed  simulations  using  a  meso¬ 
scale  model  with  simplified  physics,  and  compared  the  results  to  those  from  a  GWD 
parametrization  based  on  dry  mountain-wave  theory.  In  these  simulations,  the  relative 
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Figure  8.  The  terrain  field  (m)  for  the  domain  near  Boulder,  Colorado,  USA,  used  for  simulating  mountain  waves 
with  COAMPS®  (Coupled  Ocean/ Atmosphere  Mesoscale  Prediction  System  of  the  Naval  Research  Laboratory). 
The  three  solid  lines  denote  the  lower  boundaries  of  the  vertical  cross-sections  for  Cases  A  (06  UTC  28  December 
1998),  B  (00  UTC  31  October  1998)  and  C  (12  UTC  01  January  1999).  Abcissa  and  ordinate  are  degrees  of 

longitude  and  latitude,  respectively. 


humidity  was  set  to  1%  to  avoid  cloud  formation  and  precipitation,  the  surface  heat- 
exchange  coefficient  was  set  to  zero  to  suppress  surface  energy  transfer,  and  radiative 
forcing  was  not  allowed  in  order  to  perform  a  meaningful  comparison.  Similarly,  in 
this  study  we  limit  these  processes  in  COAMPS®  in  order  to  obtain  the  dry  mountain- 
wave  response,  except  for  a  test  case  to  see  the  impact  of  moisture  as  is  discussed  in 
subsection  4(f). 


(, b )  Dry  mountain-wave  simulations 

We  simulated  three  mountain-wave  events  over  Boulder,  Colorado,  USA  (see 
topographical  details  in  Fig.  8).  This  region  is  famous  for  frequent  occurrence  of  strong 
downslope  windstorms  along  the  lee  side  of  the  Rockies  (e.g.  Lilly  and  Zipser  1972).  We 
consider  the  innermost  nested-grid  domain  of  the  model,  which  spans  latitudes  38.70- 
41.26°N  and  longitudes  107. 63-103. 78°W.  The  vertical  domain  reaches  24  570  m  in 
height.  The  horizontal  resolution  is  2000  m  with  157  x  139  grid  points,  and  the  vertical 
resolution  varies  between  10  m  at  the  surface  and  500  m  at  the  top  with  60  levels. 

We  integrated  the  model  stalling  from  three  different  dates:  00  UTC  28  December 
1998,  00  UTC  30  October  1998  and  00  UTC  01  January  1999,  using  the  lateral  boundary 
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TABLE  1.  Details  of  the  three  case-studies  of  mountain  wave  simulations  with  COAMPS® 


Case  A 

Case  B 

Case  C 

Initialization  time 

00  UTC  28  December  1998 

00  UTC  30  October  1998 

00  UTC  1  January  1999 

Snapshot  time 

06  UTC  28  December  1998 

00  UTC  31  October  1998 

12  UTC  1  January  1999 

Low-level  wind  angle 

282 

165 

308 

from  north  (°) 

Low-level 

Westerly 

Southerly 

North-westerly 

wind  direction 

Ax  (km) 

428.5 

284.3 

514.2 

(km) 

284.3 

428.5 

514.2 

cr,,(m) 

625.1 

625.1 

625.1 

04 

0.73 

0.27 

0.34 

OC 

2.02 

2.02 

2.02 

OD 

1.43 

0.70 

1.01 

Lx 

0.48 

0.69 

0.52 

Lx 

0.69 

0.48 

0.52 

Fr0 

1.30 

1.02 

1.32 

No  (s-1) 

1.27  x  10“2 

1.44  x  10“2 

1.00  x  10"2 

Uo  (m  s-1) 

17.45 

10.36 

9.56 

m 

1.98 

1.95 

1.75 

E 

3.70 

2.71 

3.08 

Aff  (km) 

3000 

3000 

3000 

*(=  m/Xe ff)  (n-r1) 

6.58  x  10~7 

5.00  x  10“6 

5.83  x  10-7 

/tint  (m) 

6357.6 

5609.4 

5162.9 

/fief  (m) 

4783.7 

4835.6 

4794.5 

ho  (m) 

3063.5 

2673.5 

2870.9 

See  text  for  the  definitions  of  parameters.  The  height  values  are  measured  from  mean  sea  level,  while  the  surface 
level  is  at  2373.2  m. 


conditions  provided  by  NOGAPS  forecasts.  When  the  mountain  waves  were  fully 
developed,  we  examined  three  episodes:  Case  A  at  06  UTC  28  December  1998,  Case  B 
at  00  UTC  31  October  1998,  and  Case  C  at  12  UTC  01  January  1999,  where  the  low- 
level  average  wind  is  approximately  perpendicular,  parallel,  and  at  a  significant  angle  to 
the  approximate  ridge  axis  of  the  mountain  range,  respectively.  The  mountain  range  is 
oriented  approximately  north-south  in  the  centre  of  the  domain  (see  Table  1  for  details 
of  the  simulations). 

(i)  Case  A  simulation  ( 06  UTC  28  December  1 998 ).  The  terrain  and  typical  atmos¬ 
pheric  conditions  of  the  Boulder  area  (Fig.  8)  arc  generally  favourable  for  strong  wave 
generation.  At  06  UTC  28  December  1998,  the  low-level  wind  direction  is  westerly  at  the 
mountain  crest  level,  nearly  perpendicular  to  the  approximate  centre  line  of  the  ridge. 
The  winds  arc  stronger  over  the  highest  peak  of  the  mountain  range,  as  shown  in  the 
vertical  cross-section  of  the  wind  and  isentropes  (Fig.  9(a))  corresponding  to  line  A  in 
Fig.  8.  This  is  associated  with  strong  wave-breaking,  which  first  occurs  around  7  km 
over  the  steep  lee-side  slope  of  the  highest  peak  and  subsequently  propagates  down 
to  the  surface,  as  is  clearly  shown  by  steepened  isentropes  and  increased  zonal  winds 
(Fig.  9(a)).  This  downslope  windstorm  is  accompanied  by  weak  or  reversed  flow  in 
the  region  of  wave  breaking  itself  (shown  by  light  colours)  above  the  accelerated  flow 
region. 

(ii)  Case  B  simulation  (00  UTC  31  October  1998).  At  00  UTC  31  October  1998,  the 
direction  of  the  low-level  wind  is  largely  from  the  south  approximately  parallel  to  the 
mountain  ridge.  As  a  result,  the  low-level  flow  does  not  encounter  major  orography  in 
the  eastern  half  of  the  domain.  The  wave  activity  near  the  surface  is  weaker  in  Case  B 
(Fig.  9(b))  than  in  Case  A  (Fig.  9(a))  since  the  low-level  flow  is  weaker  due  to  the 
orographic  configuration  in  the  upstream  region  (see  the  southern  centre  region  along 


1906 


Y.-J.  KIM  and  J.  D.  DOYLE 


30 


25 


20 


15 


10 


Figure  9.  Vertical  cross-sections  showing  contours  of  potential-temperature  (K)  and  zonal  wind  (m  s-1,  grey 
shades — see  key)  simulated  by  COAMPS®  (Coupled  Ocean/Atmosphere  Mesoscale  Prediction  System  of  the 
Naval  Research  Laboratory)  along  the  lines  shown  in  Fig.  8:  (a)  Case  A  (06  UTC  28  December  1998);  (b)  Case  B 
(00  UTC  31  October  1998);  and  (c)  Case  C  (12  UTC  01  January  1999).  The  low-level  flow  is  from  left  to  right, 
which  is  from:  (a)  west  to  east,  (b)  south  to  north,  and  (c)  north-west  to  south-east.  The  zero  zonal  wind  contours 
are  denoted  by  long  dashed  lines.  Note  that  the  vertical  scale  of  (a)  is  different  from  those  of  (b)  and  (c). 


line  B  in  Fig.  8)  that  largely  blocks  the  low-level  flow  (clearly  seen  from  the  cross- 
section  of  the  wind  at  3000  m;  not  shown).  There  is,  however,  still  significant  wave 
activity  especially  downstream  of  the  major  mountain  peak  (Fig.  9(b)). 

(iii)  Case  C  simulation  (12  UTC  01  January  1999).  At  12  UTC  01  January  1999,  the 
low-level  winds  are  primarily  north-westerly;  thus  this  case  lies  between  Cases  A  and 
B  with  respect  to  the  low-level  wind  direction.  In  view  of  the  orographic  variations,  the 
vertical  cross-section  for  line  C  (Fig.  8)  passes  through  several  relatively  narrow  peaks 
(Fig.  9(c))  and  has  a  relatively  deep  downslope  region  in  the  downstream  half  of  the 
cross-section,  which  causes  a  fairly  strong  downslope  windstorm  over  the  lee  side  of 
the  highest  peak  in  the  centre.  The  zonal  winds  over  the  downslope  of  the  highest  peak 
are  stronger  than  in  the  surrounding  area  as  a  result  of  LLWB,  although  weaker  than 
Case  A.  Unlike  the  other  cases  the  waves  do  not  propagate  much  vertically  as  can  be  seen 
from  the  relatively  flat  isentropes  in  the  upper  part  of  Fig.  9(c).  There  are  two  distinct 
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regions  of  zonal  flow  reversal:  one  is  in  the  lee  of  the  highest  peak;  the  other  is  over  the 
downstream  end  of  the  cross-section  between  about  5.5  km  and  9  km.  shown  in  Fig.  9(c) 
inside  long  dashed  contours.  The  first  region  seems  to  be  a  direct  result  of  wave  activity 
at  low  levels,  whereas  the  second  is  due  to  the  turning  of  the  wind  direction  from  north¬ 
westerly  to  south-westerly  with  height  associated  with  a  front  along  the  downstream 
half  of  the  cross-section  (not  shown). 


(c)  Offline  evaluation  of  the  KA95+  scheme  using  the  explicit  simulations 

We  now  evaluate  the  KA95+  orographic  GWD  parametrization  scheme  using 
explicit  COAMPS®  simulations  basically  following  KA95.  The  evaluation  procedure 
consists  of  the  following  key  steps: 

•  The  horizontal  angles  of  simulated  low-level  winds  are  domain-averaged  over  the 
depth  of  2oh  ■ 

•  The  zonal,  meridional  and  vertical  components  of  the  simulated  winds  are 
decomposed  into  their  area  means  (overbar)  and  deviations  (prime)  from  the  means. 

•  The  horizontal  momentum  fluxes  arc  calculated  at  each  point,  projected  to  the 
averaged  low-level  wind  direction,  area-averaged,  and  vertically  interpolated  to  large- 

_ s 

scale  model  pressure  levels  ( pu'w '  ). 

•  The  KA95+  scheme  is  applied  to  the  area- averaged  and  vertically  interpolated 
simulated  variables,  also  using  the  orographic  statistics  calculated  from  the  topography 

- p 

used  in  the  simulations,  to  parametrize  the  momentum  flux  profile  {pu'w'  ). 

- p 

•  The  parametrized  momentum  flux  profile  {pu'w'  )  is  compared  with  the  explicitly 

_ s 

simulated  momentum  flux  profile  {pu'w'  ). 

•  In  addition,  below  the  reference  level,  the  total  parametrized  surface  drag  (GWD 
plus  BLD,  as  discussed  by  Webster  et  al.  (2003))  is  compared  with  the  explicitly 
simulated  surface  pressure  drag,  both  of  which  arc  linearly  interpolated  in  the  vertical 
and  connected  to  the  flux  profile  at  the  reference  level  (similar  to  Lane  et  al.  2000). 

(i)  Evaluation  of  Case  A  (06  UTC  28  December  1998).  We  first  evaluate  Case  A  in 
which  the  low-level  wind  is  oriented  perpendicular  to  the  mountain  ridge  axis,  which  is 
thus  a  good  test  case  for  a  quasi  2-D  situation  with  strong  LLWB  {Ftp  =  1 .30  >  Frc ;  see 
Table  1,  Fig.  9(a)).  The  low-level  wind,  which  is  the  average  between  the  surface  and  the 
height  of  2 07,  (between  2373  and  3623  m  above  sea  level,  or  1250  m  from  the  surface)  is 
westerly,  which  is  used  by  the  parametrization  at  the  reference  level.  Figure  10  shows  the 
domain-averaged  profiles  of  the  buoyancy  frequency  and  projected  wind,  which  arc  key 
input  parameters  to  the  scheme  (note  that  this  average  is  not  just  for  the  cross-section,  but 
for  the  entire  domain,  following  Broad  (1996)).  An  approximate  two-layered  structure 
is  evident,  which  is  typically  found  in  the  atmosphere  with  strong  resonant  waves.  The 
atmospheric  static  stability  for  Case  A  (Fig.  10(a))  first  increases  and  then  decreases 
with  height  in  the  lower  paid  of  the  domain,  then  rapidly  decreases  from  around  1 1  km, 
and  gradually  increases  further  up.  The  projected  wind  for  Case  A  (Fig.  10(b))  forms  a 
positive  shear  up  to  about  1 1  km  and  a  negative  shear  above  that.  Figure  1 1  shows  the 
minimum  Richardson  number  (Eq.  (8)),  which  describes  the  instability  of  the  layers  as 
diagnosed  by  the  parametrization.  This  graphically  shows  that  the  layer  instability  (Rim) 
depends  sensitively  on  its  critical  value  (we  use  Ric  =  1  as  in  LM97;  Ric  =  0.25  was 
used  by  KA95  and  K96). 

Figure  12(a)  first  compares  the  simulated  and  parametrized  momentum  flux  pro¬ 
files.  (The  ‘weighted’  parametrization  is  discussed  later  in  subsection  4(e).)  Note  that 
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Figure  10.  Profiles  for  Cases  A  (06  UTC  28  December  1998).  B  (00  UTC  31  October  1998)  and  C  (12  UTC  01 
January  1999)  of  the  domain-averaged:  (a)  Brunt-Vaisala  frequency,  and  (b)  horizontal  wind  projected  on  to  the 
direction  of  the  low-level  wind  average.  Simulations  profiled  are  by  COAMPS®  (Coupled  Ocean/Atmosphere 
Mesoscale  Prediction  System  of  the  Naval  Research  Laboratory). 


Minimum  Richardson  Number 

Figure  11.  Profiles  of  the  minimum  Richardson  number  for  Cases  A  (06  UTC  28  December  1998),  B  (00  UTC 
31  October  1998)  and  C  (12  UTC  01  January  1999).  See  text  for  details. 


the  parametrized  BLD  at  the  surface  (plus  GWD,  which  is  assumed  constant  with  height 
below  the  reference  level)  has  been  linearly  interpolated  to  the  parametrized  flux  profile 
due  to  GWD  at  the  reference  level  (denoted  by  the  horizontal  line),  while  the  surface 
pressure  drag,  calculated  in  the  direction  of  the  low-level  wind,  has  been  linearly 
interpolated  to  the  simulated  flux  profile.  Although  there  exist  some  differences  in  the 


EXTENSION  OF  OROGRAPHIC-DRAG  PARAMETRIZATION  SCHEME 


1909 


Figure  12.  Profiles  of  the  simulated  (open  circles),  parametrized  (closed  circles)  and  parametrized  with  weighted 
orographic  statistics  (triangles)  vertical  fluxes  of  the  horizontal  momentum  for:  (a)  Case  A  (06  UTC  28  December 
1998),  (b)  Case  B  (00  UTC  31  October  1998)  and  (c)  Case  C  (12  UTC  01  January  1999).  The  sign  of  the  momentum 
fluxes  is  reversed  following  the  convention  of  the  parametrization  scheme.  The  profiles  below  the  horizontal 
line  are  linear  fits  of  the  surface  pressure  drag  (simulated)  and  the  blocked-layer  drag  plus  gravity-wave  drag 
(parametrized)  to  the  reference-level  drag.  Note  that  the  horizontal  scale  of  (a)  is  different  from  those  of  (b)  and 

(c).  See  text  for  further  details. 


details,  the  matching  between  the  two  profiles  is  excellent  for  this  case.  The  use  of 
Ric  =  1  is  found  helpful  for  this  case  in  detecting  wave-breaking  at  upper  levels  where 
Rim  is  mostly  larger  than  0.25  (see  Fig.  1 1).  The  results  with  Ric  =  0.25  or  0.5  are  still 
in  reasonable  agreement  overall  (not  shown),  but  inferior  to  the  result  with  Ric  =  1  at 
upper  levels. 

(ii)  Evaluation  of  Case  B  ( 00  UTC  31  October  1 998 ).  We  next  evaluate  Case  B 
in  which  the  low-level  wind  direction  is  southerly  and  nearly  parallel  to  the  average 
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ridge  axis  (Fig.  9(b)).  This  case  is  thus  selected  to  represent  an  extreme  3-D  anisotropic 
situation.  The  flow  is  fairly  nonlinear  (Fro  =  1  -02  >  Frc,  see  Table  1)  due  to  the  charac¬ 
teristics  of  the  Front  Range  of  the  Rockies.  The  variation  of  orography  is  significant  both 
in  the  along-flow  (south-north)  and  cross-flow  (west-east)  directions  (Fig.  8).  The  hor¬ 
izontal  domain  averages  of  the  static  stability  and  projected  horizontal  wind  arc  similar 
to  Case  A,  and  show  a  two-layered  structure  with  weaker  winds  above  1 1  km  (Fig.  10). 

The  simulated  and  parametrized  momentum  flux  profiles  compare  well  (for 
Ric  =  1)  except  below  the  reference  level  (Fig.  12(b)),  and  are  not  greatly  sensitive  to  the 
choice  of  Ric  since  Ri  is  mostly  smaller  than  0.25  when  the  layer  is  unstable  (Fig.  1 1). 
It  is  noted  that  the  magnitude  of  the  flux  for  Case  B  at  the  reference  level  is  about 
an  order  of  magnitude  smaller  than  that  for  Case  A;  this  is  a  result  of  weaker  LLWB 
due  to  the  3-D  anisotropic  configuration  of  orography.  The  value  of  OD  (0.7  for  Case  B 
versus  1.43  for  Case  A,  see  Table  1)  successfully  controls  the  magnitude  of  parametrized 
GWD,  which  would  have  been  overestimated  if  OD  were  not  used.  There  is,  however, 
a  significant  difference  between  the  two  profiles  below  the  reference  level  due  to  the 
BLD  parametrization,  although  they  arc  within  the  same  order  of  magnitude  (Fig.  12(b)). 
The  parametrized  BLD  is  about  half  of  the  calculated  surface  pressure  drag. 

(iii)  Evaluation  of  Case  C  (12  UTC  01  January  1999).  We  now  evaluate  Case  C  in 
which  the  low-level  wind  is  oriented  half-way  between  Cases  A  and  B.  This  is  a  good 
test  case  for  a  situation  in  which  the  flow  is  at  an  angle  with  respect  to  the  mountain 
ridge  axis.  The  averaged  low-level  wind  is  north-westerly  and  the  flow  is  also  strongly 
nonlinear  (Fro  =  1.32  >  Frc,  see  Table  1  and  Fig.  9(c))  due  to  the  characteristics  of  this 
area  over  Rockies  (Fig.  8).  The  layer  structure  for  Case  C  is  more  complex  than  other 
cases  in  that  there  are  multiple  vertical  shear  layers:  i.e.  positive  shear  near  the  surface 
and  between  9  km  and  19  km,  and  negative  shear  elsewhere  (Fig.  10(b)).  The  overall 
structure  of  the  static  stability  is  similar  to  other  cases  (Fig.  10(a)). 

The  simulated  and  parametrized  momentum  flux  profiles  (for  Ric  =1)  show  a 
relatively  close  match.  Rim  is  very  large  at  upper  levels  (Fig.  11),  which  means 
wave-breaking  is  absent  at  upper  levels  as  confirmed  in  the  simulation  (Fig.  9(c)). 
The  magnitude  of  the  surface  flux  is  much  smaller  than  Case  A  and  similar  to  (and 
smaller  than)  Case  B,  due  to  weaker  LLWB  associated  with  orographic  anisotropy. 
The  parametrization  realistically  estimates  the  magnitude  of  the  reference-level  drag 
and  its  vertical  divergence;  however,  below  the  reference  level  BLD  is  somewhat 
underestimated  as  in  Case  B.  We  note  here  that  it  is  fairly  easy  to  retune  BLD  by  a 
factor  of  two  for  Cases  B  and  C.  The  easiest  way  is  perhaps  to  eliminate  the  factor  1  /2 
from  Eq.  (12)  or  increase  the  value  of  Cd  to  2.0  to  double  the  magnitude.  This  will, 
however,  overestimate  the  BLD  for  Case  A  by  a  factor  of  two.  We  always  try  to  prevent 
overestimation  and  allow  for  underestimation.  On  the  other  hand,  our  GWD  formulation 
shows  less  of  this  dilemma  of  overestimation  versus  underestimation,  as  seen  in  Fig.  12 
(and  later  in  Fig.  13). 

The  results  of  Cases  B  and  C  imply  that  BLD  parametrized  by  our  formulation  is 
relatively  less  dependent  on  orographic  anisotropy  and  is  somewhat  harder  to  parame¬ 
trize,  while  low-level  GWD  is  more  dependent  on  orographic  anisotropy  and  is  relatively 
easier  to  parametrize.  This  is  consistent  with  our  experience  that  an  implementation  of 
large  isotropic  BLD  into  NOGAPS  tends  to  improve  the  forecast  skill. 

(iv)  Evaluation  of  all  time  episodes  including  Cases  A,  B  and  C.  In  addition  to 
the  detailed  evaluation  of  the  above  case-studies,  the  magnitudes  of  the  reference-level 
drag  parametrized  by  the  KA95+  scheme  (Eq.  (6))  is  also  evaluated  by  simple  linear 
regression  for  all  times  taken  from  the  explicit  simulations.  These  cases  include  not  only 
wave-breaking  and  trapping  cases,  but  also  fairly  linear  cases  with  very  weak  or  virtually 
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Figure  13.  The  parametrized  gravity-wave  drag  (GWD)  at  the  reference  level  (abscissa)  versus  the  explicitly 
simulated  domain-averaged  GWD  at  the  reference  level  (ordinate).  A  total  of  twenty-six  snapshots  have  been 
taken  at  various  times  with  3-  or  6-hour  interval  from  the  mesoscale  simulations  performed  with  COAMPS®, 
including  Cases  A  (06  UTC  28  December  1998),  B  (00  UTC  31  October  1998)  and  C  (12  UTC  01  January  1999) 
as  marked  in  the  diagram.  R  is  the  linear  correlation  coefficient  of  the  least-square  curve  fit. 


no  wave  activity  (not  shown).  As  shown  in  Fig.  13,  which  consists  of  26  ‘snapshots’ 
taken  from  the  simulations,  the  correlation  between  the  parametrization  and  simulation, 
including  Cases  A,  B  and  C,  is  very  high  (with  correlation  coefficient  ~0.9). 

(d)  Flow  regime  diagram 

There  arc  different  sources  of  drag  due  to  subgrid-scale  orography  currently  con¬ 
sidered  together  with  GWD  parametrization:  BLD  due  to  low-level  flow  blocking, 
GWD  due  to  breaking  at  upper  levels  of  vertically  propagating  hydrostatic  waves, 
GWD  due  to  LLWB,  GWD  due  to  non-hydrostatically  trapped  waves  or  drag  due  to 
internal  reflection,  critical-level  absoiption,  and  orographic  lift,  etc.  Figure  1  shows 
some  of  these  sources  schematically  (see  also  Fig.  4  of  Kim  et  al.  (2003)  and  Fig.  1 
of  Zadra  et  al.  (2003)).  The  parametrizations  by  MW96  and  GSM98,  which  follow 
Shutts  (1990),  include  several  sources  of  drag  (GWD  due  to  hydrostatic  waves, 
LLWB,  trapped  lee  waves,  and  internal  wave  reflection),  while  SM00  includes  BLD, 
non-hydrostatic  and  hydrostatic  GWD;  LM97  includes  BLD  and  hydrostatic  GWD. 
The  KA95+  parametrization  includes  all  of  these,  but  GWD  due  to  hydrostatic  and 
non-hydrostatic  waves  is  treated  collectively.  The  relationship  or  harmony  among  these 
mechanisms,  at  low  levels  in  particular,  can  be  understood  using  a  ‘regime  diagram’  as 
shown,  for  example,  in  Fig.  1  of  SM00  and  Fig.  2  of  LM97.  In  a  regime  diagram,  the 
surface  pressure  drag  is  non-dimensionalized  by  its  linear  counterpart,  and  plotted  as  a 
function  of  the  inverse  Froude  number  in  order  to  identify  the  flow  regimes,  e.g.  MJ92; 
Stein  1992,  hereafter  S92;  OB96/97;  LM97;  SM00).  This  diagram  serves  as  an  effective 
measure  to  check  the  overall  behaviour  of  drag  parametrization  schemes  for  a  variety  of 
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Figure  14.  Flow  regime  diagram  that  shows  the  ratio  of  parametrized  GWD  at  the  reference  level  and/or  BLD 
to  the  linear  drag  as  a  function  of  the  domain-averaged  low-level  inverse  Froude  number,  defined  by  Eq.  (5),  but 
divided  by  OD.  Total  of  twenty-six  'snapshots’  have  been  taken  at  various  times  with  3-  or  6-hour  intervals  from 
the  mesoscale  simulations  performed  with  COAMPS®,  including  Cases  A  (06  UTC  28  December  1998),  B  (00 
UTC  31  October  1998)  and  C  (12  UTC  01  January  1999)  as  marked  in  the  diagram.  The  curves  for  MJ92,  S92  and 
LM97  have  been  reconstructed  from  Fig.  2  of  LM97  while  those  for  OB97  and  SM00  have  been  reconstructed 
from  Fig.  7(a)  of  SM00.  The  thick  solid,  long-dashed  and  short-dashed  curves  (corresponding  respectively  to 
Cases  A,  B  and  C)  are  obtained  from  the  analytic  version  of  the  parametrization  (GWD+BLD)  by  analytically 
varying  wind  and  stability,  but  with  basic  orographic  statistics  fixed  (Table  1).  See  text  for  further  details. 


flow  states.  SM00  categorized  the  flow  over  orography  depending  on  the  low-level  value 
of  the  Froude  number:  when  Fro  is  smaller  than  its  critical  value  (typically,  Frc  ~  0.8) 
the  regime  is  governed  by  linear  dynamics;  when  Fro  is  between  Frc  and  a  typical 
threshold  value  (of  less  than  2)  the  regime  is  governed  by  nonlinear  dynamics  and  the 
drag  becomes  at  least  twice  the  linear  value;  when  Frc  is  greater  than  the  threshold 
value  the  regime  is  considered  to  be  governed  by  near-surface-flow  dynamics,  while  it 
is  recognized  that  there  can  be  significant  overlap  between  the  last  two  regimes. 

Figure  14  shows  such  a  diagram  obtained  from  the  KA95+  scheme  using  all  times 
of  the  explicit  simulations,  including  Cases  A,  B  and  C,  in  comparison  with  results  from 
selected  other  studies.  Note  that  Fro  defined  by  Eq.  (5)  has  been  divided  by  OD  (Eq. 
(4))  to  make  it  comparable  to  other  studies.  For  the  linear  drag,  we  use  the  expression 
for  the  2-D  hydrostatic  non-rotating  frictionless  GWD,  KpoNoUocr^  with  k  =  8  x  10~6 
(McFarlane  1987).  Note  that  the  magnitudes  of  our  non-dimensional  drag  depend  on  our 
choice  of  k,  and  thus  the  comparison  made  with  other  studies  can  only  be  relative.  Our 
results  show  that  the  total  drag  is  clearly  divided  into  two  groups;  the  upper  one  with 
values  greater  than  2.0,  and  the  lower  one  smaller  than  2.0.  These  groups  correspond 
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directly  to  values  of  OD  which  are  greater  and  smaller,  respectively,  than  its  average 
value  of  1.08  with  four  directions.  Some  key  features  of  the  upper  group  are  also  found 
in  other  studies,  such  as  the  drastic  increase  of  the  reference-level  drag  near  the  critical 
value  of  Fro  (=  0.8)  and  gradual  decrease  of  the  drag  as  Fro  increases  further  (e.g.  LM97 
for  Co  =  2  and  S92  on  the  diagram).  In  these  results  that  involve  smaller  values  of 
Fro,  GWD  is  several  times  greater  than  BLD  as  is  also  shown  by  SM00  (see  their 
Fig.  7(a));  whereas  in  the  lower  group  results  with  smaller  non-dimensional  drag  and 
larger  values  of  Fro,  GWD  is  roughly  comparable  to  BLD.  This  shows  that  GWD  is 
dominant  for  smaller  Fro,  but  BLD  becomes  equally  dominant  for  larger  Fro •  Webster 
et  al.  (2003)  introduced  a  simple  parametrization  in  which  the  magnitude  of  BLD  is 
inversely  proportional  to  that  of  GWD,  which  is  consistent  with  our  results  except  that 
their  BLD  is  larger  than  GWD.  We  also  constructed  the  diagram  with  respect  to  the 
original  definition  of  Fro  by  Eq.  (5)  and  found  qualitatively  similar  results,  except  that 
the  lower  group  is  displaced  to  left  on  the  diagram  by  about  0.5  of  a  Froude  number 
unit  (not  shown).  Moreover,  we  found  virtually  no  correlation  between  Fro  and  OD. 
This  suggests  that  the  magnitude  of  the  drag  depends  largely  on  orographic  anisotropy, 
regardless  of  Fro-  Since  OD  enters  into  our  scheme  only  through  Fro,  the  magnitude 
of  drag  is  well  controlled  by  OD  in  the  scheme,  as  shown  in  the  comparison  with 
the  explicit  simulations  (see  Fig.  13).  Our  results  also  imply  that  there  can  be  large 
differences  between  3-D  and  2-D  orography  in  the  magnitudes  of  drag  as  shown  by 
earlier  studies,  and  also  illustrated  by  the  two  distinct  groups  of  points  in  Fig.  14. 
For  example,  MJ92,  OB 97,  FM97  (Cd  =  1)  and  SM00  (large  Fro  portion)  correspond  to 
our  3-D  results  whereas  S92,  FM97  (Cd  =  2)  and  SM00  (small  Fro  portion)  correspond 
to  our  2-D  results. 

Furthermore,  for  a  more  instructive  comparison  with  other  studies,  we  fixed  basic 
orographic  statistics  and  analytically  varied  the  wind  and  stability  for  the  three  Cases 
(A,  B  and  C,  see  Table  1),  although  a  regime  diagram  constructed  this  way  with  our 
scheme  may  not  be  very  realistic  since  some  orographic  statistics  cannot  be  analytically 
obtained.  The  analytic  drag  curves  (denoted  by  thick  curves  in  Fig.  14)  are,  overall, 
comparable  to  our  actual  simulation  results  except  for  Case  B  (thick  long-dashed  curve), 
which  is  only  about  the  half  of  the  simulated  value  (marked  by  B  on  the  diagram).  This  is 
probably  due  to  the  fact  that,  although  the  orientation  of  the  main  mountain  ridges  gives 
a  rather  small  value  of  OD  (i.e.  the  horizontal  mountain  aspect  ratio),  the  complex  nature 
of  the  orography  in  the  domain  contributes  to  fairly  nonlinear  and  vertically  propagating 
behaviour  of  the  waves  as  can  be  seen  in  Fig.  9(b),  which  is  not  well  represented  by 
the  crude  analytic  formula.  In  addition,  the  analytic  drag  curve  for  Case  A  (thick  solid 
curve)  roughly  follows  the  upper  drag  group  of  other  studies  (somewhat  smaller  between 
1.2  <  Fro  <  2.5)  except  for  Fro  smaller  than  ~0.8,  which  is  significantly  larger  than 
other  studies.  On  the  other  hand,  the  curve  for  Case  C  (thick  short-dashed  curve)  runs 
through  the  lower  group  while  that  for  Case  B  (thick  long-dashed  curve)  generally  goes 
below  the  lower  group.  It  is  evident  that  our  peak  value  of  the  analytic  drag  is  located  at 
a  smaller  value  of  Fro/OD  than  other  studies.  This  may  be  due  partially  to  the  effect  of 
OD,  and  also  to  the  mathematical  nature  of  our  formulation  that  gives  large  values  for 
small  values  of  the  static  stability  (see  Eq.  (6)). 

(e)  Weighted  averages  of  the  orographic  statistics 

In  the  introduction,  we  discussed  that  column  (i.e.  local)  GWD  parametrizations 
become  increasingly  hard  to  justify  as  model  resolution  increases  and  we  introduce  a 
non-local  version  of  the  scheme  in  appendix  A.  For  computational  reasons  under  the 
parallel  computing  architecture,  however,  our  GWD  parametrization  currently  remains 
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local.  The  sensitivity  of  OA  and  L x  to  the  location  of  orography  relative  to  the  grid  is 
in  practice  effective  for  distinguishing  different  typical  near- surface  flow  states,  but  it  is 
undesirable  in  principle:  it  would  be  better  if  the  total  drag  were  to  be  independent  of 
the  relative  location  of  orography  and  grid.  A  very  crude  way  to  alleviate  this  deficiency 
is  to  apply  global  smoothing  to  the  global  drag  field  at  the  end  of  each  time  step,  which 
may  in  some  sense  be  similar  to  using  the  non-local  version  of  the  scheme  presented 
in  appendix  A.  A  more  refined  method,  however,  is  to  take  a  weighted  average  of  the 
orographic  statistics  also  including  orography  of  nearby  grid  boxes  depending  on  the 
wind  direction.  Moreover,  the  choice  of  eight  predetermined  directions  to  calculate  OA 
and  Lx  was  mainly  motivated  by  its  simplicity,  and  also  our  experience  that  too  detailed 
wind  directions  can  excessively  constrain  the  subgrid-scale  wind  directions  that  change 
locally.  This  choice,  however,  may  invoke  a  sudden  jump  in  the  value  of  the  orographic 
statistics  with  a  small  change  in  the  wind  direction.  A  weighted  average  of  the  statistics 
will  reduce  this  effect.  Appendix  B  introduces  an  example  of  such  weighting  procedure. 

The  weighting  tends  to  smooth  the  orographic  statistics  and  slightly  increase  or 
decrease  the  aspect  ratio  depending  on  the  region.  We  compared  the  new  parametrization 
results  using  the  ‘weighted’  statistics  of  CM,  Lx  (and  hc)  with  the  original  formulations 
(see  Fig.  12).  Overall,  the  GWD  parametrization  is  similar  to  the  original  one,  while 
the  BLD  parametrization  produces  larger  differences  than  the  GWD  counterpart.  The 
BLD  parametrization  is  more  sensitive  to  the  weighting  (through  Lx)  than  the  GWD 
parametrization,  and  thus  may  require  further  calibration  for  improvement. 


(/)  Effects  of  moisture 

Durran  and  Klemp  (1982a)  investigated  the  effects  of  moisture  in  the  development 
of  trapped  mountain  waves,  and  reported  that  the  wave  response  can  be  amplified  or 
damped  depending  on  the  environment  and  the  height.  We  investigate  this  effect  by 
performing  an  additional  simulation  of  Case  A  with  moisture  included  (Fig.  15(a))  and 
comparing  it  with  the  dry  counterpart  (Fig.  12(a)).  We  find  some  noticeable  differences 
due  to  moisture  in  the  development  of  waves  such  as:  larger  wave  amplitudes,  more 
lateral  spreading  and  breaking  of  wave  packets  especially  at  upper  levels,  and  some 
differences  shown  in  the  domain-averaged  momentum  flux  profile  mostly  in  the  middle 
levels  (Fig.  15(b)).  The  parametrization  appears  to  be  less  responsive  to  vertical  drag 
divergence  between  15  and  7  km  when  moisture  is  present. 

Durran  and  Klemp  (1982b)  derived  a  modified  Brunt-Vaisala  frequency  that  in¬ 
cludes  the  moisture  contribution  to  account  for  the  effect  of  moisture.  Surgi  (1989) 
also  formulated  a  moist  Brunt-Vaisala  frequency  (Nm)  for  use  in  an  orographic  GWD 
parametrization  scheme  as  an  effective  way  of  vertically  redistributing  GWD  due  to 
moisture.  We  implemented  this  formulation  in  which  the  dry  Brunt-Vaisala  frequency 
(Aii)  is  multiplied  by  a  factor  as  follows: 

Ain  =  (1  —  e)1/2Ain  (14) 

where  e  is  0.1,  0.2  or  0.3  for  a  high-,  middle-  or  low-cloud  layer,  respectively  (Surgi 
1989).  We  found  the  use  of  Eq.  (14)  in  the  scheme  slightly  reduces/increases  the 
magnitude  of  the  drag  at  lower/upper  levels,  which  corrects  the  slightly  overestimated 
flux  at  low  levels  (Fig.  15(b)).  Our  test  result,  while  limited,  suggests  that  the  pragmatic 
assumption  made  in  virtually  all  of  the  GWD  parametrization  schemes  that  ignore  the 
effects  of  moisture  may  be  valid  only  to  a  first  order  of  approximation.  Consequently  the 
parametrizations  may  need  to  be  re-evaluated  to  fully  incorporate  the  effect  of  moisture. 
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Momentum  Flux  [N/m2] 

Figure  15.  (a)  As  in  Fig.  9(a)  for  Case  A  (06  UTC  28  December  1998),  but  with  moisture  included;  (b)  as  in 

Fig.  12(a)  (dry),  but  also  with  moisture  included,  together  with  the  moist  Brunt-Vaisala  frequency  adjusted  as 

given  by  Eq.  (14). 
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5.  Concluding  remarks 

Evaluation  of  an  orographic  drag  parametrization  scheme  in  a  large-scale  model  is 
challenging,  since  its  success  often  depends  on  its  harmony  with  other  model  physics  as 
well  as  its  own  physical  adequacy.  Evaluations  with  mesoscale  models  and  observations 
of  proven  quality  are,  therefore,  a  logical  course  of  research.  Perhaps,  a  community¬ 
wide  joint  intercomparison  effort  with  well-established  mesoscale  models  will  be  useful 
to  further  advance  the  parametrization. 

The  KA95  and  KA95+  orographic-drag  parametrization  scheme  is  basically  a 
statistical  fit  of  physically  constrained,  semi-empirical  formulations  to  mountain-wave 
simulations.  The  results  presented  in  this  study  arc  highly  encouraging.  There  remain, 
however,  issues  in  the  design  and  evaluation  of  the  extended  scheme.  As  described  in 
subsections  2(a),  4(e)  and  appendix  A,  further  study  is  required  to  fully  include  the 
non-local  effect  of  drag  under  the  parallel  computing  architecture.  Further,  as  discussed 
in  subsection  4(e),  the  sensitivity  of  the  orographic  statistics  to  the  relative  location  of 
orography  is  not  desirable  in  principle,  and  an  approximate  method  to  alleviate  this 
deficiency  has  been  presented  in  appendix  B .  The  scheme  also  needs  to  be  evaluated 
offline  for  other  regions  with  more  diverse  values  of  orographic  anisotropy  than  the  four 
values  used  in  this  study. 

In  this  study  we  did  not  investigate  the  simulated  momentum  flux  profile  below  the 
mountain  crest.  Rather,  we  compared  explicitly  calculated  surface  pressure  drag  with 
parametrized  BLD  below  the  crest  as  an  effort  to  evaluate  the  parametrization.  However, 
Laprise  and  Peltier  (1989)  discussed  as  an  extension  of  the  Eliassen-Palm  theorem  that, 
below  the  crest  of  a  mountain,  the  sum  of  the  vertical  profiles  of  Reynolds  stress — which 
is  calculated  excluding  underground  points — and  (wave)  pressure  drag  is  approximately 
constant  in  height  for  stationary  state  conditions  in  the  absence  of  lateral  boundary  fluxes 
(see  their  Fig.  19).  It  may  be  worthwhile  to  systematically  verify  this  argument  against 
explicit  simulations. 

It  has  been  justified  that  BLD  is  a  more  physical  representation  of  the  orographic 
blocking  effect  in  the  subgrid-scale  (e.g.  LM97),  compared  with  an  earlier  method  of  en¬ 
hancing  resolved-scale  orography,  such  as  envelope  or  silhouette  orography.  Orographic 
roughness  parametrizations  are  widely  used  to  represent  enhanced  surface  friction  (e.g. 
MW96;  GSM98)  although  they  may  suffer  from  deficiencies  in  some  other  applica¬ 
tions  due  to  its  temporally  invariant  nature  as  pointed  out  by  Boer  and  Lazare  (1988). 
Despite  the  inherent  differences  discussed  in  the  introduction,  orographic  roughness 
parametrizations  arc  sometimes  expected  to  perform  a  function  similar  to  BLD  (as  noted 
from  the  fact  that  it  is  often  called  form  drag).  It  is  also  found  that  some  BLD  schemes 
are  being  used  together  with  form-drag  formulations.  The  interaction  among  all  of  the 
drag  mechanisms  parametrized  near  the  surface,  such  as  friction  drag,  form  drag,  BLD, 
GWD,  drag  associated  with  convectively  generated  gravity  waves  (e.g.  Chun  et  al.  2004) 
and  even  convective  momentum  transport  (Gregory  et  al.  1997),  etc.  is  complex  and  not 
yet  clearly  understood;  accordingly  they  arc  not  fully  represented  in  large-scale  models. 
As  implied  by  the  discussion  in  the  literature  related  to  the  interpretation  of  the  regime 
diagram  as  well  as  the  relationship  between  BLD  and  form  drag  (e.g.  LM97,  GSM98, 
SM00,  Webster  et  al.  2003),  more  work  is  needed  on  including  the  effects  of  gravity- 
wave  processes  in  large-scale  models. 

The  original  KA95  scheme  is  currently  paid  of  the  NCEP  (National  Centers  for 
Environmental  Prediction)  global  forecast  models  (Alpert  et  al.  1996;  Kanamitsu  et  al. 
2002).  The  extended  scheme  (KA95+)  presented  in  this  study  is  being  tested  with 
the  ALPHA  version  of  NOGAPS  (Eckermann  et  al.  2004),  which  now  uses  the  mean 
orography  (Kim  and  Hogan  2004). 
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Figure  A.l.  A  staggered  grid  for  weighting  of  the  effective  orographic  length  Lx  (Fig.  7).  Shown  are  a  grid  box 
(domain)  centred  at  point  22  and  an  example  of  a  box  centred  at  point  1 1 ,  which  is  shifted  from  22  by  half  grid 
interval  in  both  directions.  See  text  for  further  details. 


Appendix  A 

Inclusion  of  the  non-local  effect  in  the  KA95/KA95+  GWD  parametrization  scheme 

Due  to  the  non-local  nature  of  gravity  waves  that  can  propagate  out  of  or  into 
a  grid  box  column,  the  inclusion  of  the  influence  of  the  waves  on  neighbouring  grid 
columns  may  be  an  important  issue  in  GWD  parametrizations  (see  subsection  2(a)).  For 
example,  if  the  model  grid  were  placed  as  shown  by  the  dashed  line  in  Fig.  A.  1  the  GWD 
generated  by  the  mountain  in  the  left-side  grid  domain  would  reach  the  right-side  grid 
domain  where  no  mountain  is  present.  Kim  (1992)  incorporated  this  effect  by  regarding 
the  orographic  asymmetry,  OA,  as  a  free  parameter  and  obtaining  the  expected  value  of 
the  parametrized  drag  (tgwd)  for  the  entire  range  of  OA.  The  integration  of  tgwd  over 
OA  includes  the  contribution  of  OA  due  to  orography  not  only  in  the  corresponding  grid 
cell,  but  also  in  the  neighbouring  grid  cell  along  the  direction  of  the  low-level  wind. 
Mathematically,  it  can  be  expressed  as: 

^  r  OAmax 

E(r gwd)  =  — - — —  /  tgwd(GA)  d(OA),  (A.l) 

DAmax  —  UAm in  J0Amin 

where  OA  (Eq.  (1))  is  normalized,  i.e.  OAmax  =  1  and  OAm jn  =  —  1. 

The  expected  value  of  tgwd  (Eq.  (6))  is  obtained  numerically  using  Eq.  (A.l)  for 
a  discrete  set  of  OAs  between  OA  max  and  ft4mm  (e.g.  varying  from  —1  to  +1  with 
an  interval  of  0.1)  and  also  using  the  averages  of  the  flow  parameters  (i.e.  N,  U,  and 
Ftp)  extending  over  to  the  upstream  and/or  downstream  grid  cell,  possibly  depending 
on  the  location  of  the  bulk  orography  (i.e.  the  sign  of  the  local  value  of  OA).  In  the 
example  shown  in  Fig.  A.l,  where  the  orography  is  near  the  downstream  cell  boundary 
(OA  ~  —1),  the  averages  of  the  flow  parameters  are  taken  over  the  larger  grid  cell 
extended  downstream.  A  series  of  comparisons  were  made  between  the  original  (Eq. 
(6))  and  non-local  (Eq.  (A.l))  versions  off-line,  and  also  on-line  with  a  non-parallel ized 
version  of  a  general  circulation  model  (Kim  1992).  It  was  found  that  two  versions 
produce  quantitatively  different  but  qualitatively  similar  results,  particularly  in  respect 
of  the  impact  on  the  large-scale  simulations.  Implementation  of  this  expression  into 
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a  massively  parallelized  model  code,  however,  causes  significant  degradation  of  the 
computational  efficiency  due  to  the  required  communication  with  neighbouring  grid 
cells  at  each  time  step.  Improved  optimization  methods  are  needed  to  overcome  this 
technical  difficulty. 


Appendix  B 

Weighted  averages  of  the  orographic  statistics 
In  order  to  make  the  effective  orographic  length  (LY)  less  local  (see  subsection 
4(e)),  we  consider  a  staggered  grid  shown  in  Fig.  A.l.  Lx  can  be  calculated  with 
weighting  for  the  representative  wind  directions  (Fig.  6)  as: 

W  :  <X>i(Lv)23  +  tOl(Lx)22  +  m(Lx)2 1 

S:  co\(Lx)n  +  <j)2(Lx)22  +  <*>i(Lx)32 
SW :  o>i(L.v)ii  +  a>2(Lx)22  +  &>i(Ly)33 
NW  :  CO\  (F.v)  13  +  <J>2(Lx)22  +  <yi(^x)3t, 

where  a>\  and  a>2  are  the  weighting  coefficients  (2o>\  +  coi  =  I )  for  which  we  choose 
0.25  and  0.5,  respectively.  Furthermore  hc,  which  is  calculated  similarly  to  Lx  with 
respect  to  predetermined  wind  directions  (Fig.  6),  can  also  be  weighted  following 
Eq.(B.l). 

For  taking  a  weighted  average  of  OA,  we  use  either  Eq.  (B.l)  or  the  following 
expression  based  on  Fig.  2: 

W  :  coiOAmv  +  co20A\v  +  coiOAsw 
S :  coiOAsw  +  M20As  +  coiOAse 
SW  :  a>iOA\y  +  ol>20Asw  +  co\OA$ 

NW  :  (v\OAn  +  o)20Anw  +  co\OAw, 

where  a>\  and  a>2  are  0.125  and  0.75,  respectively.  Note  that  the  weighting  coefficients 
for  OA  are  different  from  those  of  Lx,  due  to  the  difference  between  OA  and  Lx  in  the 
subgrid  boxes  considered  for  counting  the  orographic  height  points.  We  also  tested  the 
method  used  above  for  Lx  to  calculate  OA  and  found  qualitatively  similar  results.  Use 
of  these  formulations  will  make  the  orographic  statistics  more  robust. 
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